Oscillations of a Bose-Einstein condensate rotating in a harmonic plus quartic trap 
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We study the normal modes of a two-dimensional rotating Bose-Einstein condensate confined in 
a quadratic plus quartic trap. Ifydrodynamic theory and sum rules are used to derive analytical 
predictions for the collective frequencies in the limit of high angular velocities, f2, where the vortex 
lattice produced by the rotation exhibits an annular structure. We predict a class of excitations with 
frequency in the rotating frame, irrespective of the mode multipolarity m, as well as a class 

of low energy modes with frequency proportional to |m|/f2. The predictions are in good agreement 
with results of numerical simulations based on the 2D Gross-Pitaevskii equation. The same analysis 
is also carried out at even higher angular velocities, where the system enters the giant vortex regime. 

PACS numbers: 03.75.Kk, 03.75.Lm, 67.40.Vs, 32.80.Lg 



I 

O 

o 



> 
o 



o 

-I— > 



I 

o 
o 



X 
S3 



The availabiUty of traps with stronger than harmonic 
confinement opens up new scenarios in rotating uhra- 
cold gases. In principle such traps permit the realiza- 
tion of configurations rotating with arbitrarily high an- 
gular velocities, since the confining potential is always 
stronger than the repulsive centrifugal term. The first 
experiments in this direction are reported in Ref. P|. 

The stationary configurations of rotating Bose- 
Einstein condensates in the presence of a harmonic plus 
quartic trap have already been the subject of several the- 
oretical papers 0, These calculations predict novel 
vortex structures reflecting the interplay between the 
centrifugal and confining forces. In particular, if the an- 
gular velocity, fJ, exceeds a critical value, the centrifugal 
force overcomes the harmonic confinement giving rise to 
a hole in the center of the condensate. For large angu- 
lar velocities the radius of the resulting annulus increases 
linearly with fi, while the width of the annulus decreases 
like l/ri. For such geometries the dynamical behavior 
of the gas exhibits new features, whose investigation is 
the main purpose of this work. In particular, along with 
excitations involving radial deformations of the density, 
one expects the occurrence of low frequency sound waves 
propagating around the annulus. 

In this work we will calculate the frequencies of the low- 
est modes by developing an analytical description using 
hydrodynamic theory and sum rules, as well as carrying 
out simulations based upon the numerical solution of the 
Gross-Pitaevskii equation. For simplicity we will restrict 
our discussion to 2D configurations, valid for fast rotating 
condensates strongly confined in the axial direction. 

The expression for the trapping potential is given by 
the sum of quadratic and quartic components 
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^Jh/MLd±^ is the characteristic harmonic oscillator length 
where M is the atomic mass, r 



r 



is the two- 



dimensional radial coordinate and A is the dimension- 
less parameter characterizing the strength of the quartic 
term. In the following, we use dimensionless harmonic 
oscillator units, where wj. and d_L are the units of fre- 
quency and length respectively. 

When the angular velocity $7 is sufficiently high, a lat- 
tice of quantized vortices is formed. If A > 0, then for 
exceeding a critical value, 17^ > 1, the equilibrium con- 
figuration in the rotating frame corresponds to a vortex 
lattice with a hole in the center 0, Q • At even larger an- 
gular velocities the system is expected to undergo a tran- 
sition to a giant vortex where all the vorticity is confined 
to the center of the annular condensate. In the follow- 
ing we will mainly restrict the discussion to the former 
regime which is more accessible experimentally, although 
we briefly discuss the giant vortex at the end. 

For a vortex lattice the dynamics can be described by 
introducing the concept of diffused vorticity within a hy- 
drodynamic picture. This approach has already success- 
fully described the dynamics of rotating configurations 
in harmonic traps [^. Such an approximation is valid 
provided that the Thomas-Fermi condition ^ <C d is sat- 
isfied, where ^ is the healing length and d is the width 
of the annulus . In addition, the healing length should 
be small compared to the distance I — I/a/TJ between 
vortices, I. 

In the rotating frame the linearized rotational hydro- 
dynamic equations take the form 

^■dn + V -{nodv) = , (2) 



dt 
d_ 

di 



Sv + gV Sn + 2fl /\Sv = , 
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Here 



is the harmonic oscillator frequency, dj_ ~ 



where no is the equilibrium density, g is the coupling con- 
stant, and Sn and Sv are the density and velocity varia- 
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tions respectively. For an effectively 2D system, uniform 
in the axial direction over a length the coupling con- 
stant can be written a,s g — AirNa/Z where N is the 
number of particles and a is the 3D s-wave scattering 
length. The integrated density is normalized to unity. 

For > the equilibrium density in the presence of 
the potential is given by Q 



no 
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where i?i_2 are the inner and the outer radius of the an- 
nulus, respectively. The mean square radius is hence 
(r2) = / r^nodr = (i?f + i?i)/2. It is useful to in- 
troduce the variable ( = (r^ - R\/2) / where 
j_ = R2 ± Rf. Hence C varies from — 1 to 1 and is zero 
at the mean square radius of the cloud. We also recall 
that Rl = {Q^ - 1)/A and R'i = (O^ - 1)/A, where the 
angular velocity for the formation of the hole, related to 
the healing length 5, is given by il/j = (1 + 2^/\/^Y/'^ = 
+ (12A2gA)i/3 0. 

For large angular velocities, R\ 
increases quadratically with Q. while i?^ (proportional to 
the area) remains constant. Hence the radius of the an- 
nulus R+ 1 increases linearly with $7 whereas the width 
of the annulus, d = R2 — Ri, decreases like 1/Q. 

The hydrodynamic equations (O and © can be solved 
by expressing the radial and azimuthal components of the 
velocity field Sv in terms of Sn, and looking for solutions 
of the form 6n = (5n(C)e""'^e~*'^*, where m is the az- 
imuthal quantum number, is the azimuthal angle and 
uj is the excitation frequency in the rotating frame. For 

> il/i, the equation for the density becomes 
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Sn + 2mnXR^_C5n+ 
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Eq. (jS)) can be significantly simplified in the large angular 
velocity limit fi^ ^ 1 where R\ ^ f7^/A, by neglecting 
the terms of order of R^_/R\ cx l/il^. This case leads to 
a class of solutions with w cx il, obeying the equation 







(6) 



and having the form of Legendre polynomials Pj {Q , with 
j = 1, 2, . . . 'sj. The corresponding eigenfrequencies are 



= [4 + j(j + l)]f^' 



(7) 



yielding, for the most relevant j = I mode, the prediction 
UJ = V6f7. Remarkably, result ((TJ is independent of both 
the oscillator frequency uj± and the strength A of the 
quartic potential. Furthermore, it is independent of the 
value and the sign of m. The linear dependence of uj on 
n can be simply understood using the macroscopic result 



UJ = cq for the sound wave dispersion. The sound velocity 
is given by the dilute gas expression Mc^ = gn with gn cx 
Ai?i independent of ft while g cx l/d cx 
Recalling that i?^ ~ J7^/A one immediately finds w oc SI. 

Result {71 has been derived in the large limit. Solu- 
tions of Eq. (jSJ holding for all O > $7^, can be found for 
A ^ 0, where Vlh ~ 1 and the terms in R'^_/R\ cx A^/^ 
are negligible. For the j = 1 mode we find the result 
cj^ — 6ri^ — 2. When < 1 the solutions for A — > tend 
to those obtained in Ref. 4] by solving the problem with 
a rotating harmonic potential. 

The collective oscillations can also be investigated us- 
ing a more microscopic approach based on sum rules. Let 
us introduce the p-energy weighted moments 
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relative to a generic excitation operator F = J2k=i /(^)fc' 
where Eno is the energy difference between the excited 
state |n) and the ground state |0). 

A useful estimate of the frequency of the monopole 
compression mode (Ad), excited by the operator /(r) = 
r^, can be obtained using the ratio between the en- 
ergy weighted (mi) and inverse energy weighted (to-i) 
moments. The former can be expressed in terms of 
commutators as mi(F) = {[F,[H, F]]) /2 = 27V (r^), 
where H = i/kin + Hext + -ffint — flLz is the many- 
body Hamiltonian in the rotating frame with interac- 
tion term Hint — 9j2i<j^i''^i ~ '''j)- contrast, the 
inverse energy weighted moment can be calculated in 
terms of the monopole static polarizability to be r7i_i = 
— {N/M)d{r'^)/duj\ (in dimensional units), where the 
derivative should be calculated at constant angular mo- 
mentum. In the Thomas-Fermi approximation one finds 
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For f7 < ri?i, i?+ is the Thomas- Fermi radius, which can 
be found by solving the cubic equation R\{4:\R\ — 'iVl^ + 
3) = 12g/7r. For D. > Qh, since Rl = (i?f + i?f) = 
{fl^ — 1)/A, one finds the simple result uj = \/6f7^ — 2, 
which is consistent with the hydrodynamic prediction for 

The result of estimate as a function of 17, is re- 
ported in Fig.Q] We compare to the numerical results ob- 
tained by solving the 2D time dependent Gross-Pitaevskii 
equation, where the numerical methods are detailed in 
Ref. Starting from the stationary solution, the mode 
is excited by a sudden change in the confining poten- 
tial, which, after some short time, is reset to its original 
form. The subsequent changes in the radius are then 
analyzed to extract the frequencies of oscillation. We 
have performed simulations at g = 1000 for A = 0.5 and 
A = 10~^. The latter value of A is similar to that used in 
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FIG. 1: Frequency of the lowest m = mode as a function of 
the angular velocity Q (in units of the trap frequency) for g — 
1000. The sum rule estimates Q for A = 0.5 and A = 10"^ 
are plotted as solid black and dotted lines respectively, while 
the results of solving the GP equation numerically are plotted 
as solid and open circles. The gray line is the asymptotic 
prediction uj — and the dashed lines represent the critical 
frequencies for hole formation, ^Ih, for both values of A. 



the experiments of Ref. 1], where the numerical solution 
of the linearized hydrodynamic equations for the lowest 
monopole oscillation at fl < flh is also presented. Fig. 
shows that the sum rule approach provides an excellent 
estimate of the monopole frequency. Fig. also reveals 
a cusp in the mode frequency at = flh- This behavior 
results from the use of the Thomas- Fermi approximation 
for in Eq. (jnj and is smoothed out by quantum effects 
included in the full GP solution. 

Sum rules can also be applied to excitations of the 
form / — T-I^le*™"^, which carry multipolarities different 
from zero. We consider the moments mf{F) = mi{F) + 
mi{F^) and m^{F) = m3{F) + m3{F^). These can be 
easily expressed in terms of commutators as mf{F) = 
{[F\[H,F]]) andm+{F) = {[[F^,H],[H, [H,F]]]). From 
the previous expressions for the dipole (m — 1) and 
quadrupole (m — 2) operators V and Q, and using the 
Thomas-Fermi approximation 7] for > flh one pre- 
dicts the following results for the ratio between the cubic 
and energy weighted sum rules 



m+(X>) 
m+(Q) 



= 517^ - 1 + - 
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In the large limit Eqs. ifTUI) and |[TT|) both yield y^Q, 
for the excitation frequency, which does not coincide with 
the prediction of Eq. (UJ. This result, which is inconsis- 
tent with a one-mode assumption, reveals the existence 
of additional modes not described by Eq. 0. Indeed, as- 



suming that is exhausted by the j = 1 modes, the fact 
that the ratio m J / m]*" is smaller than the corresponding 
frequency ^/6il implies that mf includes contributions 
from lower frequency modes. In particular one concludes 
that the latter account for 1/6 of the total rrii moment. 

The V, dependence of these lowest frequency modes can 
be simply inferred from Eq. jSJl where, neglecting higher 
order corrections, one finds that the frequency should 
be proportional to \m\XRt/^- These modes can be in- 
terpreted as describing a sound wave directed along the 
azimuthal direction, in contrast to the high-lying modes 
which correspond to a radial shape oscillation of the an- 
nulus. The coefficient of proportionality can be estimated 
from the ratio between the energy and inverse energy 
weighted sum rules. The low-lying modes contribute only 
1/6 of the moment; the m~^^ sum rule, which is ex- 
pected to be exhausted by the low lying modes, is given 
by the static response x- In the large 17 limit, the lin- 
earized hydrodynamic equations with a multipole pertur- 
bation give mti = -X = {Ntt/ g)R'i{n'^ /2Xy"^i Since 
~ 2iVm2(f72/2A)l"'l-i in the same limit, we find the 
frequency uj = {mf /6mti)'^/^ = {V2/6)\m\XRl/n. The 
same result can also be derived using a variational anal- 
ysis of Eq. (O. It is also worth noticing that lj cx A^/'^ 
tends to zero in the A ^ limit. 

Fig. 12 shows a comparison between the analytical and 
numerical results for the high-lying and low-lying dipole 
and quadrupole modes. One sees good agreement be- 
tween the two datasets at high f2, validating the sum 
rule approach used here. In the numerical simulations, 
the high-lying modes depart from the ^/E^l dependence 
for Q < il/i. The behavior at small f2 is qualitatively 
similar to the one exhibited in a rotating harmonic trap, 
where only one mode per branch is present. In particu- 
lar for A ^ 1 and < 1 the equations of rotational h; 
drodynamics in the rotating frame give the result ^ 
tj(m = ±2) = -\/2 — 17^ ip Q for the two quadrupole fre- 
quencies, while for the dipole one has ui{m = ±1) = l=Fi7. 
At large fi, the numerical results also show that the low- 
lying quadrupole mode frequency is larger than that of 
the dipole mode by a factor of two, in agreement with 
the arguments presented above. 

The excitation energies in the laboratory frame are re- 
lated to those in the rotating frame by -Eiab — -E'rot +mfl. 
For a proper identification of the modes in the lab frame, 
it is crucial to consider the sign of the azimuthal quantum 
number m associated with each excitation. For this pur- 
pose it is useful to evaluate the strengths cr+ = |(n|i^|0)p 
and a~ = |(n|F'l'|0)p, relative to the operators F and i^^ 
exciting states with angular momentum ±m. A careful 
analysis of the response function reveals that the upper 
quadrupole level corresponds to an m = —2 mode, the 
m = +2 strength relative to this level being extremely 
small. A different situation takes place for the low-lying 
level. When 17 < 1 this level has mainly an to = +2 char- 
acter, as in the case of Ref. 0. For il > il/j, instead, both 
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FIG. 2: Numerical results for the high- lying (circles) and low- 
lying (squares) \m\ — 1 (open) and |m| — 2 (closed) mode 
frequencies as a function of the angular velocity Q (in units 
of the trap frequency), for g = 1000 and A — 0.5. These are 
compared to the analytical results, where solid lines show the 
sum rule estimates for the low-lying modes while the gray line 
is the asymptotic prediction uju — \/^Q, expected to hold at 
large fl. The vertical dashed line denotes the value of Qh- 

the m = ±2 strengths significantly differ from zero. For 
example, at f2 = 5 and for the parameters of Fig. |21 the 
numerical simulation shows that cTjt ~ 0, cr^ — 12.07V, 
— 19.97V, — 7.8N, where cr^ ^ are the strengths 
associated with the high and low-lying m — ±2 modes. 
In conclusion, we predict that in the lab frame for Q > Qh 
one should observe two m = —2 modes with frequencies 
Wfj ^ 257 and jwL — 2f2|, and one m — +2 mode with fre- 
quency ti>L + 2r2, where wh,l are the high and low- lying 
frequencies in the rotating frame (21] . We also notice that 
at high angular velocity the m = — 2 mode with frequency 
[wl — 2Q,\ is energetically unstable in the lab frame. Sim- 
ilar results are found for the dipole modes. 

We finally discuss the case of the giant vortex equi- 
librium configuration, where the velocity field of the con- 
densate is irrotational. In this case, linearizing the Gross- 
Pitaevskii equation in the rotating frame gives two cou- 
pled equations for the density and the phase variations 
6n and SS 

|j„+(!5iL-n)^ + v^(„„vjs) = o, (12) 

where Vi„ — v jr for a giant vortex with circulation v [^. 
From these equations one can derive an equation similar 
to Eq. lO) , but for the phase rather than the density. For 
large f2, the solutions are again Legendre polynomials, 
but with eigenfrequencies — + 1)0^ where j > 1. 
Hence the j ~ I mode has the same frequency for both 
the irrotational and solid body cases, but the frequencies 



for the j > 1 modes are different. In the case of the 
low-lying modes for m ^ 0, using the sum rule or hydro- 
dynamic methods discussed earlier, we find a frequency 
that has the same l/fi dependence as in the vortex lattice 
case, but is larger by a factor 3^/^. 

In summary, we have studied normal modes of a Bose 
condensate in a harmonic plus quartic potential using an- 
alytic methods (hydrodynamic equations and sum rule) 
and numerical solution of the Gross-Pitaevskii equation. 
At large angular velocities fl we find a radial mode with 
a frequency ^/6fl independent of the mode multipolarity 
and value of A, as well as low-lying modes corresponding 
to waves around the annular condensate. 

This research was supported by the Ministero 
dellTstruzione, deU'Universita e della Ricerca. Some of 
the work was performed at the Kavli Institute for Theo- 
retical Physics, supported in part by the National Science 
Foundation under Grant No. PHY99-0794. 



[1] V. Bretin, S. Stock, Y. Seurin, and J. Dalibard, 
Phys. Rev. Lett. 92, 050403 (2004); S. Stock, V. Bretin, 

F. Chevy, and J. Dahbard, Europhys. Lett. 65, 594 (2004). 
[2] A.L. Fetter, Phys. Rev. A 64, 063608 (2001); U.R. Fis- 
cher and G. Baym, Phys. Rev. Lett. 90, 140402 (2003); 

G. M. Kavoulakis and G. Baym, New J. Phys. 5, 51.1 
(2003); E. Lundh, Phys. Rev. A 65, 043604 (2002); 
A.D. Jackson, G.M. Kavoulakis, and E. Lundh, Phys. Rev. 
A 69, 053619 (2004); A.D. Jackson and G.M. Kavoulakis, 
Phys. Rev. A 70, 023601 (2004). 

[3] A.L. Fetter, B. Jackson, and S. Stringari, 

cond-mat/0407119 
[4] M. Cozzini and S. Stringari, Phys. Rev. A 67, 041602(R) 

(2003). 

[5] The j — solution is rejected due to its unphysical nature. 
Indeed, for the case m — this solution does not conserve 
the particle number, while for m 7^ it gives rise to a 
divergent velocity field at the condensate boundaries. 

[6] The resuh lj'^ =60.^-2 can also be found within the 
hydrodynamic formalism by solving Eq. ^ perturbatively 
with respect to l/fi. 

[7] For the dipole operator one finds mf{T>) — 2N and 
m^{V) = 2N{l + 3n'^ +4X{r'^)), while for the quadrupole 
one instead has mf{Q) = SN{r'^) and m^{Q) = 
167V[(6!^^ -I- l){r2) + + 3A(r^) - 6fl{£.)], where = 
-{Oydx'^ + d^/dy^) and = -id/d(p. In the Thomas- 
Fermi diffused vorticity approach the equilibrium velocity 
is vo ^ ft Ar, so that (p^) = n^{r^) and (4) = ^(r^). 

[8] The result for the quadrupole frequencies has been con- 
firmed in the experiment of P.C. Haljan, I. Coddington, 
P. Engels, and E.A. Cornell, Phys. Rev. Lett. 87, 210403 
(2001), where a harmonic trap was employed. 

[9] This scenario is based on the assumption that in the ro- 
tating frame only two levels (high and low-lying) are ex- 
cited by the operators F and . Numerical integration of 
Eq. actually shows the occurrence of a very small split- 
ting between the low-lying m = ±2 modes, which however 
barely affects the main results of the present analysis. 



